Transcriptional analysis in multiple barley varieties identifies signatures of waterlogging response

Abstract Waterlogging leads to major crop losses globally, particularly for waterlogging‐sensitive crops such as barley. Waterlogging reduces oxygen availability and results in additional stresses, leading to the activation of hypoxia and stress response pathways that promote plant survival. Although certain barley varieties have been shown to be more tolerant to waterlogging than others and some tolerance‐related quantitative trait loci have been identified, the molecular mechanisms underlying this trait are mostly unknown. Transcriptomics approaches can provide very valuable information for our understanding of waterlogging tolerance. Here, we surveyed 21 barley varieties for the differential transcriptional activation of conserved hypoxia‐response genes under waterlogging and selected five varieties with different levels of induction of core hypoxia‐response genes. We further characterized their phenotypic response to waterlogging in terms of shoot and root traits. RNA sequencing to evaluate the genome‐wide transcriptional responses to waterlogging of these selected varieties led to the identification of a set of 98 waterlogging‐response genes common to the different datasets. Many of these genes are orthologs of the so‐called “core hypoxia response genes,” thus highlighting the conservation of plant responses to waterlogging. Hierarchical clustering analysis also identified groups of genes with intrinsic differential expression between varieties prior to waterlogging stress. These genes could constitute interesting candidates to study “predisposition” to waterlogging tolerance or sensitivity in barley.

different aspects of plant growth, but at the same time, they contribute to the onset of a multifaceted response. One of the most important changes caused by waterlogging is the reduction in oxygen availability (hypoxia) to roots and also to shoots in the case of flooding (Bailey-Serres, Fukao, et al., 2012;Bailey-Serres, Lee, & Brinton, 2012;Sasidharan et al., 2017). Other parameters that change and may affect plants negatively upon waterlogging and flooding include the release of toxic chemical compounds in the soil (e.g., iron) (Setter & Waters, 2003), reduced availability of important nutrients such as nitrates (Setter & Waters, 2003), changes in the soil microbiome (Hartman & Tringe, 2019;Wang et al., 2017), or a decrease in light availability because of flood water turbidity.
In recent years, the study of waterlogging stress and of the resulting hypoxic stress has led to the discovery of essential and conserved oxygen-sensing mechanisms in plants (reviewed in Hammarlund et al., 2020;Holdsworth & Gibbs, 2020).
An important oxygen-sensing pathway in plants requires a set of oxygen-dependent plant cysteine oxidase (PCO) enzymes that oxidize the N-terminus of proteins starting with a cysteine residue (Weits et al., 2014;White et al., 2018White et al., , 2017, including a set of group VII ethylene response factor (ERF-VII) transcription factors (White et al., 2018(White et al., , 2017) that act as master regulators of the hypoxia response program (reviewed in Giuntoli & Perata, 2018). Under normal oxygen conditions, these ERF-VII transcription factors are targeted for degradation by the ubiquitin-dependent N-degron pathway, following oxidation of their N-terminal cysteine residue by PCO enzymes (Gibbs et al., 2011;Licausi et al., 2011;Weits et al., 2014). In contrast, under hypoxic conditions, the activity of PCOs is inhibited, thus preventing the oxidation of ERF-VIIs' N-terminal cysteine and their subsequent N-degron-mediated degradation. As a result, under hypoxic conditions, ERF-VIIs accumulate in the nucleus, where they bind to cis-regulatory elements to activate the hypoxia response program. Notably, the promoters of conserved hypoxia response gene families are not only enriched for cis-regulatory elements bound by transcription factors of the ERF-VII family (e.g., the HRPE motif; Gasch et al., 2016) but also for motifs bound by basic helix-loop-helix (bHLH), MYB, and WRKY transcription factors (Reynoso et al., 2019), thus suggesting the involvement of other families of transcription factors in the regulation of the hypoxia response program (Lee & Bailey-Serres, 2021). This is in agreement with the upregulation of transcription factor-coding genes belonging to the different families mentioned above as part of hypoxia response in plants (Mustroph et al., 2010). Cross-comparisons of transcriptomic datasets from different plant species further allowed the identification of core genes of the hypoxia response program (Mustroph et al., 2010(Mustroph et al., , 2009Reynoso et al., 2019). These include genes involved in (i) the regulation of carbon metabolism in order to facilitate adenosine triphosphate (ATP) production via glycolysis and nicotinamide adenine dinucleotide (NAD+) regeneration through fermentation pathways and (ii) the regulation of important signaling pathways (e.g., mitogen-activated protein kinase [MAPK] signaling) and molecules (e.g., reactive oxygen species [ROS]) to promote plant tolerance and survival upon hypoxia.
Numerous transcriptomic analyses using either flooding or waterlogging treatments have been conducted with the model plant Arabidopsis thaliana (see list in Brazel & Graciet, 2023). These studies revealed that core aspects of the transcriptional reprogramming upon waterlogging or flooding are shared with the response to hypoxia van Veen et al., 2016). Furthermore, the transcriptomic comparison of eight different natural accessions of Arabidopsis, which had been previously identified as being either sensitive or tolerant to flooding , indicated that the transcriptional response of roots and shoots differ (van Veen et al., 2016). This work identified sets of shoot-or root-expressed genes that might contribute to the flooding-tolerance phenotype of some natural Arabidopsis accessions (van Veen et al., 2016), which could be relevant to improve crop tolerance to waterlogging/flooding. Notably, the transcriptional reprogramming in response to hypoxia and flooding is also accompanied by other genome-wide regulatory mechanisms such as epigenetic changes, chromatin remodeling (Reynoso et al., 2019), and changes in translation (Lee & Bailey-Serres, 2019;Reynoso et al., 2019).
Waterlogging is an important source of crop losses; however, not all crops are equally affected (de San Celedonio et al., 2014;Kaur et al., 2020). Barley is particularly sensitive to waterlogging with reported crop losses of up to 20-25% (Byrne et al., 2022;de San Celedonio et al., 2016Liu et al., 2020;Miricescu et al., 2021;Setter & Waters, 2003). These losses are particularly severe if waterlogging occurs at the heading stage (de San Celedonio et al., 2016Liu et al., 2020;Setter & Waters, 2003) but also at early developmental stages (de San Celedonio et al., 2016. The multifaceted nature of waterlogging stress, as well as the complexity of plant responses to this stress, have made it difficult to identify marker genes of waterlogging tolerance, as well as target genes that could be modified to improve crop tolerance to waterlogging. Nevertheless, in recent years, genetic approaches, as well as linkage mapping and genome-wide association studies, have led to the identification of potential quantitative trait loci (QTLs) and target genes to improve barley tolerance to waterlogging (Bertholdsson et al., 2015;Borrego-Benjumea et al., 2020;Broughton et al., 2015;Gill et al., 2019;Li et al., 2008;Manik et al., 2022;Zhang et al., 2017Zhang et al., , 2016Zhou, 2011).
In agreement with the complex response to waterlogging, these QTLs have been identified based on a wide range of phenotypic traits, including leaf yellowing, chlorophyll fluorescence, yield traits, adventitious root formation, aerenchyma formation, and ROS levels. QTLs relating to root traits, including the ability to form adventitious roots and aerenchyma under waterlogged conditions (Manik et al., 2022;Zhang et al., 2016), may be of particular relevance considering their known roles in facilitating water and nutrient uptake as well as in oxygen distribution during waterlogging (Zhang et al., 2015).
Despite the recent progress made in the identification of QTLs for waterlogging tolerance and the realization that transcriptomics may be used to identify potential candidate genes relevant to waterlogging tolerance Reynoso et al., 2019;van Veen et al., 2016), only few such studies have been conducted with barley. Recent studies have characterized the transcriptomic response of waterlogging-tolerant and waterlogging-sensitive varieties (Borrego-Benjumea et al., 2020;Luan et al., 2022). From these datasets, the authors identified a handful of genes that could be of importance for waterlogging-tolerance mechanisms in barley (Luan et al., 2022). In another study, a comparative proteomic analysis of one waterlogging-sensitive variety and one waterlogging-tolerant variety identified proteins that accumulate differently in waterlogging-sensitive or waterlogging-tolerant germplasms .
Here, to dissect the transcriptomic response of barley to waterlogging, we selected two 2-row and two 6-row winter barley varieties based on the differential expression of known hypoxia response marker genes after waterlogging treatment. These varieties were chosen from the Association Genetics Of UK Elite Barley (AGOUEB) population (Thomas et al., 2014) and from the list of recommended barley varieties in Ireland (where the study was conducted; list at the time the experiments were carried out). Our transcriptomic study focused on roots because of their central role in mediating waterlogging tolerance (Zhang et al., 2015). We also included the model spring barley variety Golden Promise to provide a reference dataset to the community (this variety is widely used to target specific genes using molecular approaches). We identified sets of common genes that are consistently differentially expressed in barley in response to waterlogging. Hierarchical clustering identified groups of genes with intrinsic differential expression between varieties prior to waterlogging stress.
Low or elevated expression of some of these genes could "predispose" some varieties to waterlogging tolerance or sensitivity. In sum, the datasets presented serve as an additional reference for the study of waterlogging response in barley and provide insights into potential avenues of research to improve waterlogging tolerance in this crop.

| Plant material
Cultivars used in this study (Table S1) included winter varieties that originated from the AGOUEB population of barley (Thomas et al., 2014), as well as the model spring variety Golden Promise (obtained from Teagasc, Oak Park, Ireland) and a winter variety Infinity (obtained from Teagasc, Oak Park, Ireland), which was on the recommended list of barley varieties in Ireland (where the study was conducted) at the time of the initial field trial whose results were taken into account for varietal selection (Byrne et al., 2022). Selected varieties for further characterization, including transcriptome profiling, were Golden Promise (spring variety, two-row), Regina and Infinity (winter varieties, two-row), and Passport and Pilastro (winter varieties, six-row).

| Plant growth conditions
Plants were grown in a plant growth room under long-day conditions (16 h light/8 h dark) at 15 C (constant temperature), approximately 45% relative humidity. Light intensity was determined to be $138 μmol/m 2 /s and was provided by LED bulbs (Philipps LED tubes High Output, T8 20W/865).

| Soil preparation and seed germination
Commercial John Innes No. 2 (Vitax, UK) soil was soaked in water after filling round pots of 9 cm diameter and 9 cm height without compressing the soil. Untreated seeds were sown directly in soil at a depth of 2 cm. The sown seeds were stratified in the dark for 14 days at 4 C to ensure homogenous germination. After cold treatment, the pots were transferred to the plant room for germination and growth.

| Waterlogging
Plants were grown as indicated above for 10 days (corresponding to L1/L2 stage) for transcriptomic experiments and for 10-14 days (corresponding to L1/L2 stage) for phenotypic characterization. The pots were then transferred to a large tub, which was subsequently filled with tap water up to 1 cm above soil level. The water level was kept constant throughout the duration of the experiment (Figure 1a). The plants were kept in the same growth conditions and were treated for 15 days for phenotypic characterization or for shorter periods of time, as indicated in the text. Control plants were left in the same growth conditions but received normal watering (every 2 days, avoiding any standing water in the trays). For the recovery period, plants were taken out of the water and kept in the same growth conditions with a normal watering regime.

| Height measurements
Plant height was taken from the soil surface to the tip of the tallest leaf at the indicated time points.

| Total RNA extraction
Roots of plants grown under control conditions or subjected to waterlogging were rinsed under running tap water and frozen in liquid nitrogen. The tissue was ground in liquid nitrogen, and total RNA was extracted using Spectrum™ Plant Total RNA kit (Sigma-Aldrich). For each condition (waterlogged or untreated), the root systems of three plants of the same variety were pooled prior to grinding for total RNA extraction. This experiment was conducted independently at least three times to obtain samples from at least three biological replicates, as indicated in the figure legends.

| Reverse transcription coupled to quantitative polymerase chain reaction (RT-qPCR)
Total RNA was reverse transcribed using an oligo (dT) 18 primer and the RevertAid reverse transcriptase (Thermo) according to the manufacturer's instructions. For reverse transcription reactions, 1 μg of total RNA was used. cDNA obtained was used for qPCR with a Lightcycler 480 (Roche). Each PCR reaction mix contained 5 μL of 2Â SYBR green master 1 (Roche), 1 μL cDNA, 1 μL of 10 μM primers, and 3 μL of molecular biology grade water. LightCycler melting curves were obtained to check for single peak melting curves for all amplification products. The second derivative maximum method was used to analyze the amplification data. The resulting Cp values were converted into relative expression values using the comparative Ct method (Livak & Schmittgen, 2001). All primer sequences are provided in Table S2. One reference gene (HvACTIN) was used to normalize the RT-qPCR data following the screening of a set of reference genes during waterlogging ( Figure S1).

| Next-generation sequencing of RNA samples
For RNA-seq analysis, waterlogging treatment was applied as outlined above for 24 h. For each condition (waterlogged or untreated), the root systems of three plants of the same variety were pooled prior to grinding for total RNA extraction. This experiment was conducted independently three times to obtain samples from three biological replicates (i.e., for each variety, six RNA samples were sent for sequencing, corresponding to three biological replicates for the untreated plants and three biological replicates for the waterlogged plants). Following RNA extraction, RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent). All RNA samples had RNA integrity (RIN) values >7.0. Library preparation and single-ended 50 bp next-generation sequencing was performed by BGI Genomics (Hong Kong) using the DNBseq sequencing platform.  Raw RNA-sequencing reads were aligned to Morex V3 using bowtie2 (v2.4.5) (Langmead & Salzberg, 2012). Files were converted from .sam to .bam files and indexed using samtools (v1.15.1) (Danecek et al., 2021). Gene abundance was estimated using stringtie (v2.1.7) (Pertea et al., 2015) (Table S3). Differential gene expression analysis was performed using the Bioconductor package DeSeq2 (Love et al., 2014) in R (Team, 2022) using a design in which the factors Variety and Treatment were combined into a single factor to model multiple condition effects. The results of multiple comparisons were extracted and filtered by adjusted p-value < 0.05 (Table S4). Principal component analysis (PCA) was performed using pcaExplorer (v2.20.2) (Marini & Binder, 2019) in RStudio (v2022.02.0 + 443). Gene ontology (GO) analysis was performed using ShinyGO (v0.75c) (Ge et al., 2020).

| Analysis of RNA-seq data
Read count values were transformed by variance stabilizing transformations to normalize for library size and composition (Table S5). Venn diagrams were generated using InteractiVenn (Heberle et al., 2015).
The means of variance stabilizing transformations read counts from three biological replicates generated from DeSeq2 were filtered for k-means clustering as follows. Deseq2 differentially expressed gene (DEG) analysis was performed, and results for the 25 comparisons shown in Figure S6A were extracted. A list of all 11,613 DEGs filtered by adjusted p-value < 0.05 was generated by combining DEGs from waterlogged to control samples from the same variety, DEGs from each control to every other control sample, and DEGs from each waterlogged to every other waterlogged sample. Next, the DEGs with a mean of <10 normalized reads were removed leaving 10,882 genes for clustering analysis. Clustering analysis was performed using the k-means function in R (v3.6.2) with the arguments centers = 25, nstart = 1000, iter.max = 300, and algorithm = "Lloyd" (gene cluster assignment annotated in Table S4). Clustering heatmaps were generated using ComplexHeatmap (v2.10.0) (Gu et al., 2016) in R.
To compare the datasets obtained to a previously published RNA-seq dataset, raw sets were downloaded from NCBI's Gene Expression Omnibus (GSE144077) (Borrego-Benjumea et al., 2020) and Sequence Read Archive (PRJNA889532) (Luan et al., 2022). Raw RNA-sequencing reads were aligned to Morex V3 as described above.
Differential gene expression analysis of waterlogging versus control for each variety was performed on 0 h control and 24 h waterlogging datasets for the varieties Franklin and TX9425 from Luan et al. (2022) as described above. The same differential gene expression analysis was also performed on 72 h control and 72 h waterlogging datasets for the varieties Deder 2 and Yerong from Borrego-Benjumea et al. (2020).
To compare the datasets to a previously published list of 28 genes found in QTLs for waterlogging tolerance, we downloaded the gene list from tab. S3 in Manik et al. (2022). This gene list used Morex V1 (r1) gene IDs. To compare this list to our data, we performed a BLASTN search of the coding sequence (CDS) of each of these genes to the Morex V3 genome using Phytozome 13 (Goodstein et al., 2012) and used the transcript with the highest percentage identity to a Morex V3 transcript for the comparison (Table S4).

| Characterization of the transcriptional response of different barley varieties using hypoxia response marker genes
The previous characterization of 403 varieties from the AGOUEB collection under waterlogged and control conditions in the field (Byrne et al., 2022) provided some information on the differential physiological response of these varieties to waterlogging, while also highlighting the difficulties associated with the study of waterlogging tolerance under field conditions. Based on this initial study, we selected a subset of 20 varieties that behaved differently under waterlogged conditions, with the aim to assess their transcriptional response to waterlogging under controlled conditions. The variety Golden Promise was also included as a model (spring) barley variety that is widely used to generate targeted mutations through Agrobacterium-mediated transformation. After identifying HvACTIN as a suitable reference gene in waterlogged roots ( Figure S1A and S1C), we first carried out a time course experiment with four selected varieties and determined the relative expression of HEMOGLOBIN1 (HB), ALCOHOL DEHYDROGE-NASE 1 (ADH1), and PYRUVATE DECARBOXYLASE 1 (PDC1). These genes are hypoxia response markers commonly used to monitor waterlogging response (Loreti et al., 2020;Mendiondo et al., 2016;Mustroph et al., 2010). This initial analysis indicated that (i) the expression of these three hypoxia response genes peaked around 24 h after the onset of the waterlogging treatment and (ii) there were differences between the four varieties in terms of the amplitude of the transcriptional upregulation ( Figure S1B). For example, at 24 h after the beginning of the waterlogging treatment, HvHB was expressed at higher levels in Pilastro compared with Arma, Tapir, and Masquerade. In addition, the upregulation of HvHB and HvADH1 was stronger than that of HvPDC1, possibly making these first two genes more suitable to identify varieties with differential transcriptional regulation of the hypoxia response program.
We next tested the expression of HvHB, HvADH1, and HvPDC1 at 24 h of waterlogging in the set of 21 barley varieties we selected based on Byrne et al. (2022) (Figure 1a). As expected, under normal watering conditions, most varieties had low levels of expression for each of the hypoxia response genes, with the exception of Pilastro, which exhibited higher HvHB expression (Figures 1b,c and S2). Waterlogging stress triggered the upregulation of the hypoxia response genes in all varieties tested, but some differences in their response were observed. As previously, differences were more marked with HvHB and HvADH1 (Figure 1b-d) than with HvPDC1 ( Figure S2), making these first two genes more suitable to identify varieties with a differential transcriptional response to waterlogging. After 24 h of waterlogging, some varieties reached higher expression levels than the average relative expression observed for these genes in the population of 21 varieties. These included Pilastro and Regina for both HvHB and HvADH1, as well as Dura, Vesuvius, Siberia, Mahogany, and Tamaris for HvADH1. In contrast, some varieties showed reduced upregulation of hypoxia response genes compared with the population average. For example, varieties such as Golden Promise, Passport, Isa, Louise, and Retriever had lower expression levels of both HvHB and HvADH1 compared with the average of the population.
Other varieties, such as Infinity, had average expression levels of HvHB and HvADH1.
Based on these results, Pilastro and Regina were chosen as representative six-row and two-row varieties, respectively, that potentially exhibit a stronger transcriptional response to waterlogging, while Passport was selected as a six-row variety that had a more dampened transcriptional response. We also included Infinity as a two-row variety for further characterization because it is on the recommended list in Ireland (where this study was carried out) and has an average transcriptional response to waterlogging. As highlighted above, Golden Promise was included as a reference variety.

| Effect of waterlogging on the growth of selected varieties
The growth of the five selected varieties was characterized under controlled conditions following a two week waterlogging treatment and a six week recovery period with normal watering. Height measurement showed that the growth of Infinity and Passport was more affected by waterlogging than that of Pilastro, Regina, and Golden Promise, whose height was mostly unaffected by the treatment (Figure 2a). In addition, Golden Promise, Infinity, and Pilastro produced fewer tillers, whereas Regina's tiller number was largely unaffected (Figure 2b). Because root traits have been shown to be important for waterlogging tolerance (Zhang et al., 2015), we determined the length of the primary root, as well as the number of seminal roots after F I G U R E 2 Effect of waterlogging on shoot and root growth. (a) Plant height after 14 days of waterlogging and a 6 week recovery period. For Golden Promise, Passport and Pilastro, three biological replicates (four plants per biological replicate) were carried out. For Regina and Infinity, two biological replicates (four plants per replicate) and three biological replicates (four to nine plants per replicate) were carried out, respectively. (b) Number of tillers after 14 days of waterlogging and a 6 week recovery period. Three biological replicates were carried out with at least four plants per replicate. (c) Length of the primary root after 14 days of waterlogging. For Golden Promise, Passport, Pilastro, and Regina, three biological replicates with two to seven plants each were carried out. For Infinity, five biological replicates were carried out with three to four plants per replicate. (d) Number of seminal roots after 14 days of waterlogging. For Golden Promise, Passport, and Pilastro, three biological replicates (two to four plants per replicate) were carried out. For Regina, two biological replicates were conducted (four plants per replicate). For Infinity, four biological replicates were carried out (four plants per replicate). In (a)-(d), data from different biological replicates are color coded. Error bars correspond to standard deviations. Statistical analysis: multiple unpaired t-test with multiple comparison correction (Holm-Sidak method). Asterisks: *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001. ns, not statistically significant.
2 weeks of waterlogging. The length of the primary root was reduced upon waterlogging stress for all varieties except Regina (Figure 2c), while the number of seminal roots increased in all varieties and the difference was statistically significant for all except Infinity (Figure 2d).

This phenotypic characterization suggests that a variety such as
Regina, which has a stronger transcriptional response for HvADH1 and HvHB, appears to be less impacted by waterlogging than other varieties.

| Overview of transcriptional responses to waterlogging in root tissue of selected varieties
To determine the transcriptional responses of Pilastro, Regina, Passport, Infinity, and Golden Promise to waterlogging, we subjected these varieties to waterlogging for 24 h. This time point was chosen based on the time course experiment described above that showed peak expression of HvHB, HvADH1, and HvPDC1 at 24 h of waterlogging.
Total RNA was extracted from roots for RNA-seq analysis. Reads obtained were aligned to the Morex barley genome (version 3; Table S3 and Figure S3A-S3C), and differential gene expression (adjusted p-value < 0.05) was determined for each of the waterlogged varieties relative to the corresponding untreated variety (Table S4) Figure S3D); or (iii) whether they were two-row or six-row varieties (PC3: 11.34% variance; Figure S3E). In this analysis, the varieties Golden Promise, Infinity, and Passport could also be separated from Pilastro and Regina ( Figure S3F; PC4: 9.92% variance). This suggests that there are underlying differences between these two groups of varieties that are not explained by a known variable.

The number of DEGs identified in the different varieties varied
considerably: for Golden Promise and Passport, we identified 220 and 592 DEGs, respectively, while Infinity, Pilastro, and Regina had a substantially higher number of DEGs (1307, 1512, and 1233, respectively). In all varieties tested, the proportion of upregulated genes was higher than that of downregulated genes (Figure 3b). We verified in our datasets the expression of the HvHB, HvADH1, and HvPDC1 genes whose upregulation was initially used to monitor the transcriptional response of the different varieties to waterlogging by RT-qPCR ( Figure 1). Induction of all three genes upon waterlogging was found ( Figure 3c). We further identified additional homologs of HvADH1 and HvHB and confirmed that these homologs were also upregulated in response to waterlogging ( Figure S4).  waterlogging (Luan et al., 2022). In addition, Borrego-Benjumea et al.

| Genome-wide transcriptional reprogramming of barley in response to waterlogging
(2020) analyzed the transcriptomic response of waterlogged barley roots in 2 week old Yerong and Deder2 for 72 or 120 h. Both varieties had been chosen because of their tolerance to waterlogging (Li et al., 2008;Takeda, 1989). Raw RNA-seq data were analyzed using the same pipeline as for our datasets. PCA analysis of all datasets showed that most of the variation was due to the origin of the dataset ( Figure S5A; PC1: 68.39% variance), likely because waterlogging stress is known to vary considerably based on experimental differences, such as soil type. As expected, the different samples could also be separated based on treatment ( Figure S5A;  Figure S5G).

| Identification of unique expression signatures between varieties
To compare the responses to waterlogging between the five varieties used in this study, k-means clustering was performed on 10,882 genes identified as differentially expressed in (i) at least one variety when comparing waterlogged to control samples for each variety or (ii) between the control samples (i.e., these have intrinsic differential expression between at least two varieties in the absence of waterlogging) or (iii) between the waterlogged samples from different varieties (i.e., genes that are differentially expressed between at least two varieties under waterlogging) ( Figure S6A). Clustering analysis revealed distinct gene expression patterns between samples ( Figure 5a). A number of clusters contained genes that were upregulated in multiple varieties in response to waterlogging compared with untreated plants (i.e., clusters 4, 7, 8, 14, and 21). The largest cluster (#14) comprised 802 genes (Table S4) with a centroid expression that was higher in all varieties in response to waterlogging (Figures 5b and   S6). Cluster 14 included homologs of known hypoxia response genes r3.3HG0252910) ( Figure S4A and S4B). Furthermore, transcription factors from the ERF; bHLH; basic leucine zipper (bZIP); no apical meristem, ATAF1,2, and CUC2; and WRKY families were also identified in cluster 14 (Table S4)  In contrast, clusters 5, 9, 13, 15, 20, and 24 comprised genes that were downregulated in response to waterlogging in multiple varieties.
Cluster 9 contained 435 genes that were downregulated in response to waterlogging (Figures 5d and S6), including 14 of the 16 genes that were downregulated in waterlogged samples from all varieties ( Figure 4c and Table S4). Among these genes was an aminotransferase (HORVU.MOREX.r3.2HG0190680), whose top-scoring Arabidopsis homolog is γ-aminobutyric acid transaminase (AT3G22200). In Arabidopsis, the expression of this gene is also reduced in response to hypoxia (Branco-Price et al., 2008). Another member of cluster 9 whose expression was downregulated in all varieties in response to waterlogging is a cation exchanger (HORVU.MOREX. homeostasis and was shown to be downregulated in response to hypoxia and waterlogging (Wang et al., 2016). Multiple transcription factors from the heat-shock and MYB transcription factor families were also found in cluster 9. This is consistent with the observation that differentially expressed MYB transcription factors were predominantly found to be downregulated in response to waterlogging in Yerong and Deder2 (Borrego-Benjumea et al., 2020). GO analysis of the 435 genes in cluster 9 showed an enrichment for terms related to lipid and carbohydrate metabolism, transmembrane transporter activities, and cell periphery (Figure 5e), similarly to the GO terms identified in all genes downregulated in response to waterlogging (Figure 4e).
Clustering analysis also identified gene expression signatures that were specific to individual varieties: Cluster 1 (257 genes) contained genes whose expression was the highest in the barley varieties Golden Promise and Infinity and intermediate in Passport. In contrast, the expression of these genes was markedly lower in the varieties Pilastro and Regina (Figures 5f and S6). GO analysis of genes in cluster 1 identified terms related to catabolic and metabolic processes (Figure 5g), which could be relevant to waterlogging tolerance/sensitivity. In contrast, genes in cluster 16 (359 genes) had an opposite behavior, in that they were more highly expressed in Pilastro and Regina but were only expressed at low levels in the other varieties (Figures 5h and S6). GO analysis of genes in cluster 16 revealed an enrichment for genes associated with diverse processes found in other clusters or sets of genes, except the term "nicotianamine metabolic process" which was only found in the GO analysis of cluster 16 genes. This GO term relates to metal (including iron) homeostasis, which could be of relevance to waterlogging stress tolerance because of the increased uptake of metals and their toxicity to plants in waterlogged conditions.
Notably, genes in clusters 1 and 16 were not differentially expressed in response to waterlogging. Instead, it is their intrinsic (or varietal-specific) expression level that differs between the samples.  et al., 2022, 2023). This overlap, together with the results of GO and RT-qPCR analyses, validates our datasets. We also identified a set of 98 waterlogging response genes that were common to the five datasets generated in this study. Many of these common DEGs are homologs of "core hypoxia response genes" that were identified across a range of plant species (Mustroph et al., 2009) and that are involved in processes such as cell wall metabolism, ethylene biosynthesis and signaling, carbohydrate metabolism, and ROS regulation. Multiple transcription factors belonging to the ERF, bHLH, and WRKY families, which are known to be important for the onset of the waterlogging response program, are also among these 98 common DEGs. Hence, our data indicate that, as expected, the gene regulatory network that controls waterlogging response is conserved among monocots and dicots (Mustroph et al., 2010;Tamura & Bono, 2022).

| DISCUSSION
Clustering analysis of the different DEGs identified groups of genes that behave similarly between the five varieties in response to waterlogging (e.g., all genes in clusters 14 and 9 are upregulated or downregulated, respectively). This analysis also revealed genes that are differentially expressed between the varieties in the absence of treatment (see clusters 1 and 16; Figure 5f,h). While their expression does not typically change in response to waterlogging, it is possible that either the low or high intrinsic expression of some of these genes may "predispose" certain varieties to either sensitivity or tolerance to waterlogging. In other words, these genes could constitute varietyspecific waterlogging susceptibility or tolerance factors. Validating this possibility would require (i) an in-depth characterization of the physiological responses and yield losses of a large number of varieties that would be first grouped based on the intrinsic expression level of genes of interest; (ii) dissecting the function of these genes, many of which have paralogs in barley; and (iii) generating barley plants that are mutated or that constitutively express these different genes, followed by a detailed characterization of the physiological response and yield under waterlogged conditions. Potential trade-offs on other traits of agronomic relevance and on resilience to other abiotic and biotic stresses would also need to be assessed. Considering the variety of gene functions present in clusters 1 and 16, it is however difficult to predict which genes (or gene families) are most likely to be of potential relevance to crop improvement. Their association with genomewide association studies and QTL identification would most likely help to pinpoint the best candidates. For example, we compared our data to a list of 28 genes located in QTLs for waterlogging tolerance (Table S4)  r3.7HG0736590). This gene may be relevant because regulation of potassium flux during waterlogging has been shown to be important (Gill et al., 2018). Although the expression does not necessarily change in response to waterlogging, the intrinsic expression differences between varieties may be of interest to breeding waterloggingtolerant varieties in barley.

AUTHOR CONTRIBUTIONS
AM, AJB, and EG conducted experiments; AJB, AM, JB, FW, and EG analyzed RNA-seq datasets and wrote the manuscript.

ACKNOWLEDGMENTS
The authors are thankful to William Thomas for sharing seeds of the AGOUEB population and to Susanne Barth and Carl Ng for helpful discussions. We also thank Daniel Daly for technical assistance with some of the waterlogging experiments. Open access funding provided by IReL.

CONFLICT OF INTEREST STATEMENT
The authors declare that they have no conflicts of interest.

PEER REVIEW
The peer review history for this article is available in the supporting information for this article.

DATA AVAILABILITY STATEMENT
The datasets supporting the results of this article are included within the article (and its supplemental files). RNA-seq data have been deposited with the Gene Expression Omnibus (GEO) repository (at http://www.ncbi.nlm.nih.gov/) under GSE220532.